// -*- C++ -*-------x-----------------------------------------------------------
//
//
// Class:           HenikoffProfile
//
// Author:          Enrico Negri
//
// Project Name:    BiocompUP Align
//
// Date:            02/2008
//
// Description:     Calculate a frequency profile or PSSM using Henikoff
//                  weighting scheme. Some explanations can be found in:
//
//                  Guoli Wang, Roland L. Dunbrack jr.
//                  Scoring profile-to-profile sequence alignments.
//                  Institute for Cancer Research, Fox Chase Cancer Center,
//                  Philadelphia, Pennsylvania 19111, USA. March 16, 2004.
//
// -----------------x-----------------------------------------------------------

#include <HenikoffProfile.h>

namespace Biopool
{

// CONSTRUCTORS:

HenikoffProfile::HenikoffProfile() : Profile()
{ }


HenikoffProfile::HenikoffProfile(const HenikoffProfile &orig) : Profile(orig)
{
	copy(orig);
}


HenikoffProfile::~HenikoffProfile()
{ }


// OPERATORS:

HenikoffProfile&
HenikoffProfile::operator = (const HenikoffProfile &orig)
{
	if (&orig != this)
		copy(orig);
	return *this;
}


// MODIFIERS:

void
HenikoffProfile::copy(const HenikoffProfile &orig)
{
	Profile::copy(orig);

	aliWeight.clear();
	for (unsigned int i = 0; i < orig.aliWeight.size(); i++)
	{
		vector<double> tmp;
		for (unsigned int j = 0; j < orig.aliWeight[0].size(); j++)
			tmp.push_back(orig.aliWeight[i][j]);
		aliWeight.push_back(tmp);
	}
}


HenikoffProfile*
HenikoffProfile::newCopy()
{
	HenikoffProfile *tmp = new HenikoffProfile(*this);
	return tmp;
}


// HELPERS:

void
HenikoffProfile::pCalculateWeight(Alignment &ali)
{
	const string residue_indices = "ARNDCQEGHILKMFPSTWYV";


	// --------------------------------------------------
	// 1. Calculate the positions of the first aminoacids
	//    for every sequence
	// --------------------------------------------------

	vector<unsigned int> firstAminoPos;

	firstAminoPos.push_back(0); // master sequence

	for (unsigned int j = 0; j < (numSeq - 1); j++) // other sequences
		for (unsigned int i = 0; i < seqLen; i++)
			if (ali.getTemplatePos(i, j) != '-')
			{
				firstAminoPos.push_back(i);
				break;
			}


	// --------------------------------------------------
	// 2. Calculate the positions of the last aminoacids
	//    for every sequence
	// --------------------------------------------------

	vector<unsigned int> lastAminoPos;

	lastAminoPos.push_back(seqLen-1); // master sequence

	for (unsigned int j = 0; j < (numSeq - 1); j++) // other sequences
		for (unsigned int i = seqLen; i > 0; i--)
			if (ali.getTemplatePos(i-1, j) != '-')
			{
				lastAminoPos.push_back(i-1);
				break;
			}


	// --------------------------------------------------
	// 3. Calculate weights for master sequence
	// --------------------------------------------------

	vector<double> wSeq;
	vector<unsigned int> seqSubset;

	for (unsigned int i = 0; i < seqLen; i++)
		if (ali.getTargetPos(i) == '-')
			wSeq.push_back(0.00);
		else
		{
			// Calculate the subset of sequences

			seqSubset.clear();

			seqSubset.push_back(0);
			for (unsigned int k = 0; k < (numSeq - 1); k++)
				if (ali.getTemplatePos(i, k) != '-')
					seqSubset.push_back(k+1);

			unsigned int subsetSize = seqSubset.size();


			// Calculate Cleft and Cright

			unsigned int Cleft = firstAminoPos[seqSubset[0]];
			unsigned int Cright = lastAminoPos[seqSubset[0]];

			for (unsigned int k = 1; k < subsetSize; k++)
			{
				if (firstAminoPos[seqSubset[k]] > Cleft)
					Cleft = firstAminoPos[seqSubset[k]];
				if (lastAminoPos[seqSubset[k]] < Cright)
					Cright = lastAminoPos[seqSubset[k]];
			}


			double sum = 0.00;

			for (unsigned int p = Cleft; p <= Cright; p++)
			{
				// Calculate the number of different aminoacids

				unsigned int Ndiff = 0;

				for (unsigned int index = 0; index < 20; index++)
				{
					if ((seqSubset[0] == 0) && (ali.getTargetPos(p) == residue_indices[index]))
					{
						Ndiff++;
						continue;
					}
					for (unsigned int k = 1; k < subsetSize; k++)
						if (ali.getTemplatePos(p, seqSubset[k]-1) == residue_indices[index])
						{
							Ndiff++;
							break;
						}
				}


				// Calculate the number of the same aminoacid

				unsigned int n = 0;

				if (seqSubset[0] == 0)
					n++;
				for (unsigned int k = 1; k < subsetSize; k++)
					if (ali.getTemplatePos(p, seqSubset[k]-1) == ali.getTargetPos(p))
						n++;


				sum += (1 / (double)(Ndiff * n));
			}

			wSeq.push_back((1 / (double)(Cright - Cleft + 1)) * sum);
		}

	aliWeight.push_back(wSeq);


	// --------------------------------------------------
	// 4. Calculate weights for other sequences
	// --------------------------------------------------

	for (unsigned int j = 0; j < (numSeq - 1); j++)
	{
		wSeq.clear();

		for (unsigned int i = 0; i < seqLen; i++)
			if (ali.getTemplatePos(i, j) == '-')
				wSeq.push_back(0.00);
			else
			{
				// Calculate the subset of sequences

				seqSubset.clear();

				if (ali.getTargetPos(i) != '-')
					seqSubset.push_back(0);
				for (unsigned int k = 0; k < (numSeq - 1); k++)
					if (ali.getTemplatePos(i, k) != '-')
						seqSubset.push_back(k+1);

				unsigned int subsetSize = seqSubset.size();


				// Calculate Cleft and Cright

				unsigned int Cleft = firstAminoPos[seqSubset[0]];
				unsigned int Cright = lastAminoPos[seqSubset[0]];

				for (unsigned int k = 1; k < subsetSize; k++)
				{
					if (firstAminoPos[seqSubset[k]] > Cleft)
						Cleft = firstAminoPos[seqSubset[k]];
					if (lastAminoPos[seqSubset[k]] < Cright)
						Cright = lastAminoPos[seqSubset[k]];
				}


				double sum = 0.00;

				for (unsigned int p = Cleft; p <= Cright; p++)
				{
					// Calculate the number of different aminoacids

					unsigned int Ndiff = 0;

					for (unsigned int index = 0; index < 20; index++)
					{
						if ((seqSubset[0] == 0) && (ali.getTargetPos(p) == residue_indices[index]))
						{
							Ndiff++;
							continue;
						}
						for (unsigned int k = 1; k < subsetSize; k++)
							if (ali.getTemplatePos(p, seqSubset[k]-1) == residue_indices[index])
							{
								Ndiff++;
								break;
							}
					}


					// Calculate the number of the same aminoacid

					unsigned int n = 0;

					if ((seqSubset[0] == 0) && (ali.getTargetPos(p) == ali.getTemplatePos(p, j)))
						n++;
					for (unsigned int k = 1; k < subsetSize; k++)
						if (ali.getTemplatePos(p, seqSubset[k]-1) == ali.getTemplatePos(p, j))
							n++;


					sum += (1 / (double)(Ndiff * n));
				}

				wSeq.push_back((1 / (double)(Cright - Cleft + 1)) * sum);
			}

		aliWeight.push_back(wSeq);
	}




}


void
HenikoffProfile::pCalculateRawFrequency(vector<double> &freq, double &freqGap,
	Alignment &ali, unsigned int i)
{
	if (aminoAcidOneLetterTranslator(ali.getTargetPos(i)) != XXX)
		freq[aminoAcidOneLetterTranslator(ali.getTargetPos(i))] += aliWeight[0][i];
	else
		freqGap++;

	for (unsigned int j = 0; j < (numSeq - 1); j++)
		if (aminoAcidOneLetterTranslator(ali.getTemplatePos(i, j)) != XXX)
			freq[aminoAcidOneLetterTranslator(ali.getTemplatePos(i, j))] += aliWeight[j+1][i];
		else
			freqGap++;
}


void
HenikoffProfile::pConstructData(Alignment &ali)
{
	if (!gap)
	{
		ali.purgeTargetInsertions();
		seqLen = ali.getTarget().size();
	}

	pCalculateWeight(ali);

	gapFreq.reserve(seqLen);
	for (unsigned int i = 0; i < seqLen; i++)
	{
		vector<double> tmp;
		for (AminoAcidCode i = ALA; i <= TYR; i++)
			tmp.push_back(0.00);
		profAliFrequency.push_back(tmp);
		gapFreq.push_back(0.00);
	}

	profAliFrequency.reserve(seqLen);
	for (unsigned int i = 0; i < seqLen; i++)
	{
		pCalculateRawFrequency(profAliFrequency[i], gapFreq[i], ali, i);

		double frequencySum = 0.00;
		for (AminoAcidCode j = ALA; j <= TYR; j++)
			frequencySum += profAliFrequency[i][j];

		for (AminoAcidCode j = ALA; j <= TYR; j++)
			profAliFrequency[i][j] /= frequencySum;
	}





	setSeq(ali.getTarget());
}

} // namespace
